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Recent work by the authors led to the development of a mathematical theory dealing with 'second- 
order hyperbolic Fuchsian systems', as we call them. In the present paper, we adopt a physical 
standpoint and discuss the implications of this theory which provides one with a new tool to tackle 
the Einstein equations of general relativity (under certain symmetry assumptions). Specifically, we 
formulate the 'Fuchsian singular initial value problem' and apply our general analysis to the broad 
class of vacuum Gowdy spacetimes with spatial toroidal topology. Our main focus is on providing a 
detailed description of the asymptotic geometry near the initial singularity of these inhomogeneous 
cosmological spacetimes and, especially, analyzing the asymptotic behavior of timelike geodesies 
— which represent the trajectories of freely falling observers — and null geodesies. In particular, we 
numerically construct Gowdy spacetimes which contain a black hole-like region together with a flat 
Minkowski-like region. By using the Fuchsian technique, we investigate the effect of the gravitational 
interaction between these two regions and we study the unexpected behavior of geodesic trajectories 
within the intermediate part of the spacetime limited by these two regions. 



I. INTRODUCTION 



building on earlier 
led to a mathemat- 



Recent work by the authors 
and pioneering investigations 
ical theory of the so-called second-order hyperbolic 
Fuchsian systems. From a physical standpoint, sup- 
pose that we have a system of evolution equations that 
describes the dynamics of some physical configuration. 
As it is often the case in practice, one is not able to find 
exact analytical solutions of these equations and, instead, 
one seeks a description of the dynamics in certain asymp- 
totic regimes of interest. Such an effective description is 
often found by neglecting certain terms in the evolution 
equations which, according to physical intuition or other 
formal arguments, turn out to be inessential. In many 
applications, this leading-order behavior leads one to a 
singular problem. 

In such a context, the second-order hyperbolic Fuch- 
sian theory, discussed in the present paper, allows one 
to address the following issues. First of all, it gives pre- 
cise conditions on the data and the equations whether 
the leading-order description above is actually consistent 
with the evolution equations in a well-defined sense and, 
hence, whether the intuitive or heuristic understanding 
of the physical system can be validated. It allows us to 
formulate a singular initial value problem based on 
this leading-order description and, most importantly in 
the physical applications, to actually compute approxi- 
mations of arbitrary accuracy of general solutions of the 
evolution equations beyond the leading-order description. 
This singular initial value problem is analogous to the 
(standard) initial value problem classically formulated 



for nonlinear hyperbolic equations. The leading-order 
expansion plays the role of the free Cauchy data and is 
hence referred to as asymptotic data. We thus con- 
struct solutions to the equations which have a prescribed 
singular behavior. Specifically, keeping in mind our ob- 
jective to provide suitable tools for the physical applica- 
tions, we discuss two relevant approximation schemes in 
the present work which are useful for different purposes. 
On one hand, the approximation scheme introduced in 
[l| (cf. also the earlier work can be used to compute 
numerical solutions of arbitrary accuracy and, therefore, 
is useful for quantitative statements. On the other hand, 
another important approximation method can be used to 
generate formal expansions of arbitrary order and there- 
fore provides a useful tool for qualitative studies, see also 
5]- 

The proposed theory (together with some forthcoming 
generalizations, see e.g. [Ill ]) has been found to apply to 
a variety of problems arising in physics and, especially, in 
general relativity. In earlier work, we considered Gowdy 
spacetimes with spatial toroidal topology and we applied 
the theory to the construction of the so-called asymp- 
totically velocity dominated solutions [l2|,[l3| of the 
Gowdy equations [14l - ll6| . In the present paper, we con- 
tinue this analysis and seek for a deeper physical under- 
standing of the vicinity of the singularity existing in such 
spacetimes. In particular, we study the behavior of freely 
falling observers, i.e. timelike geodesies, and we demon- 
strate that the Fuchsian method allows us to construct 
families of such curves which "start" on the singularity 
at prescribed locations. Furthermore, we describe their 
leading-order behavior. 
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We emphasize that the issues discussed in the present 
paper are motivated by the ongoing and very active re- 
search on the dynamics of inhomogeneous cosmologies; 
the reader is referred to the contributions [T^-HH for fur- 
ther details. 

The paper is organized as follows. In a first part, we 
begin with a general discussion of the second-order hyper- 
bolic Fuchsian theory and our notion of the singular ini- 
tial value problem. To this end, we first recall some basic 
material about the (standard) initial value problem for 
nonlinear hyperbolic equations. The singular initial value 
problem is discussed next and the similarities with the 
(standard) initial value problem are stressed. We list the 
main conditions which need to be checked in order to val- 
idate the proposed formulation. This allows to conclude 
whether the singular initial value problem is well-posed 
and, hence, that the proposed leading-order description is 
consistent with the given equations. Then, we outline the 
description of the two approximation schemes mentioned 
earlier in this introduction. In the second part of this pa- 
per, we apply the theory to (vacuum) Gowdy spacetimes. 
We first summarize some now classical properties of such 
spacetimes, and then move on to the core of the present 
work, that is, the discussion of the asymptotic behavior 
of freely falling observers in the vicinity of the singularity. 
Since the geodesic equation is "only" a system of ordi- 
nary differential equations (ODE's), the discussion there 
highlights the essential aspects of the Fuchsian techniques 
without being distracted by the rather technical issues 
arising for partial differential equations (PDE's). Our 
discussion demonstrates how the theory can be applied, 
how the singular initial value problem works, and what 
limitations should be kept in mind in the applications. 
We complete this second part of the paper with exten- 
sive numerical experiments leading to the construction 
of Gowdy spacetimes with Cauchy horizons. The paper 
closes with a concluding section. 



II. SECOND-ORDER HYPERBOLIC FUCHSIAN 
EQUATIONS 

A. The initial value problem 

Recall that the function u(t, x) describing the displace- 
ment from an equilibrium position u(t, x) = of a (violin, 
say) string satisfies the linear wave equation 



□ i 



dt 2 



dx 2 



= 0, 



(1) 



in which c represents the speed of sound, t is the time 
variable, and the spatial variable x varies in some inter- 
val. This equation will serve as a model equation in all 
of the following discussion; recall that, in particular, Ein- 
stein's field equations imply, in certain gauges, nonlinear 
wave equations which describe the dynamics of the gravi- 
tational field [ID, HIl ■ Before we focus on such wave-type 
equations, let us, however, study some of the principles of 



the simple linear model Eq. ([1) first. The initial value 

problem associated with the wave equation is posed as 
follows. If we choose (smooth, say) free data functions 
denoted here by u*(x) and it**(x), then there exists pre- 
cisely one (smooth) solution u{t, x) of the wave equation 
(fTJ) with the property that 



it(0, x) = u*(x) and dtu(0,x) 



t (x). 



The initial value problem associated with the wave equa- 
tion is hence well-posed, as we will say. The interpre- 
tation of the well-posedness of the initial value problem 
for the wave equation is as follows. Since the state of the 
string at the initial time uniquely determines the state of 
the string at all times, the physical theory describing the 
evolution of string via the wave equation is determinis- 
tic. Indeed, the mathematical notion of well-posedness 
of the initial value problem is strongly related to the 
physical notion of determinism and is hence a concept 
of fundamental importance. 

Note that we are ignoring here the issue of the formu- 
lation of the boundary conditions and, throughout this 
paper and without further notice, simplify the discussion 
by assuming periodicity in space. Moreover, note that 
one can write the solutions of the wave equation explic- 
itly in terms of the data functions it, and it*». We will 
not make use of this since, later, we will in fact be inter- 
ested in general nonlinear equations for which no explicit 
solution formulas are known in general. 

Let us discuss the following important interpretation 
of the initial value problem for the wave equation. For 
this consider the Taylor expansion in t of the solution at 
the initial time 

u(t, x) = u(0, x) + d t u(0, x)t + i<9 t 2 u(0, x )t 2 + ... 

The first two terms are determined by the data. All 
higher-order terms, however, are determined by the wave 
equation from the initial data, as demonstrated by the 
expansion 



u(t, x) = u*(x) + u :t ^(x)t 



+ 2 C u *\ x ) t + g c "♦♦W* + 



(2) 



Hence, the initial value problem for the wave equation 
allows us to fix freely the short-time behavior of the so- 
lutions, i.e. the time for which t<l when terms of order 
t 2 etc. are negligible. However, this means that the so- 
lution is fixed for all times, and in particular we are not 
able to control the long-time behavior in addition. 

This just observed fundamental phenomenology for the 
simple wave equation carries over to a much larger class 
of equations, namely to general hyperbolic systems^, 



1 To be precise, we mean symmetric hyperbolic systems here, after 
a reduction to first-order form. 
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also referred to as nonlinear wave equations. We will not 
give a formal definition here, cf. 24]; for the purpose of 
this paper we can mostly think of equations of the form 
Eq. (TTJ) with certain nonlinearities. Particular examples 
will be given later. The main fact is that the initial value 
problem is well-posed in the same way as it is for Eq. (fTJ) . 
In particular, the short-term behavior of solutions can be 
prescribed directly by means of free data functions, while 
the equation, as soon as the data have been prescribed, 
leaves no further freedom to influence the long-time be- 
havior. 

We note that solutions of general nonlinear hyperbolic 
equations often show severe phenomena which occur af- 
ter "longer" evolution times, for example blow up of so- 
lutions, shocks, loss of uniqueness and bifurcations, etc. 
These nonlinear phenomena are extremely important for 
many applications in modern physics and mathematics. 
The questions how and under which conditions those de- 
velop from smooth initial conditions is often particularly 
challenging. 

A key result in general relativity is that Einstein's 
vacuurro field equations imply a system of (very com- 
plicated) nonlinear wave equations plus constraint equa- 
tions. The associated initial value problem is well-posed. 
However, its formulation is more involved than the one 
met with a standard system of hyperbolic equations. On 
one hand, Einstein equations are geometric equations in 
nature and, consequently, the type and character of the 
resulting partial differential equations depend crucially 
on the choice of the coordinate gauge and the chosen 
formulation. On the other hand, Einstein equations do 
not lead to a standard initial value problem due to the 
presence of constraints. Thanks to the fact that the con- 
straints propagate, the essential properties of the initial 
value problem are, however, preserved. Indeed, the well- 
posedness of the initial value problem for the Einstein 
equations was established first, in 1952, by Choquet- 
Bruhat (H. 



B. The singular initial value problem 

As already pointed out in the introduction of this pa- 
per, many physical applications give rise to an effective 
leading-order description at some initial time, say t = 0, 
which is singular in nature. It is thus desirable to seek for 
a formulation of the initial value problem when leading- 
order terms, that are more general than the truncated 
Taylor expansions Eq. @, are prescribed. One may won- 
der whether it is possible to formulate such a singular 
initial value problem, for which we are allowed to pre- 
scribe the behavior at the vicinity of the singularity, at 
least for short times, in the same way as the (standard) 



initial value problem allows us to fix the behavior close 
to the "regular" initial time. Whenever this is possible, 
such a theory gives us the opportunity, for example, to 
study "how smooth conditions arise from singularities" , 
as opposed to "how singularities develop from smooth 
conditions" as is usually done with the (standard) initial 
value problem. 

Let us, however, mention the following difficulty first. 
It is not reasonable to expect that a singular initial value 
problem as above is well-defined for general equations, 
except possibly in the physically uninteresting set-up of 
analytic data and solutions. In practice, attention are 
concentrated to equations of a particular type and, specif- 
ically, we focus here on the class of second— order hy- 
perbolic Fuchsian PDE's introduced in [lj, that is, 



u tt (t,x) 



2a(x) 



1 



u t (t,x) 



b(x) 



u(t, x) 



t ' ' ' i 2 ~ v "'~ 7 (3) 

= t~ 2 f(t,x,u,u x ,u t ) + c 2 {t,x)u xx (t,x). 



See 11[ for the generalization to quasilinear hyperbolic 
equations. For simplicity, let us suppose for all of what 
follows that the spatial domain is one-dimensional and 
that all functions under consideration are periodic in the 
spatial variable x. The vector-valued map u is the un- 
known in ([3]), while the given matrix- valued coefficients 
a, b, c are assumed to be smooth and, for the sake of sim- 
plicity in this paper, diagonal. The source-term / is the 
nonlincarity and will be required later to satisfy a certain 
"decay condition" . This latter condition will imply that 
the terms on the left-hand side impose the leading-order 
behavior in t. In short, second-order hyperbolic Fuch- 
sian PDE's are systems of hyperbolic equations contain- 
ing singular coefficients at t = 0. We are interested in 
general solutions defined for t > and in their asymp- 
totic behavior for t — > 0+. 

The following presentation will become more transpar- 
ent if we now multiply Eq. ([3]) by the factor t 2 (which 
vanishes on the singularity) and we introduce the singu- 
lar operator D = tdt, so that 



D 2 u(t, x) + 2a(x)Du(t, x) + b(x)u(t, x) 
= f(t,x,u,tu Xl Du) + t 2 c 2 (t,x)u xx (t,x). 



(4) 



Let us emphasize that D 2 u stands for td t (td t u). The 
left-hand side of the above equation is referred to as the 
Fuchsian principal part, while the right-hand side is 
referred to as the Fuchsian source-term. Let us also 
introduce the associated operator 



L[u] :=D 2 u(t, x) + 2a(x)Du(t, x) + b(x)u(t, x) 
- t 2 (?(t,x)u xx {t,x). 



(5) 



2 Similar statements can be made in the presence of matter fields. 
In this paper, however, we restrict attention to the vacuum case. 



Now, in order to formulate the singular initial value 
problem we look for solutions u of the form 

U = Uq + W, 

where the remainder w must be of "higher order" (in t 
and at t = 0) than the leading— order term uq. This 
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will be described in more details below; see also [l| for 
precise mathematical statements. Provided, for a given 
Mo, a unique remainder w exists, which is smooth for all 
t > 0, such that u is a solution, then the singular initial 
value problem associated with uq will be called well- 
posed. The function uo plays the role of the (in general) 
"singular data" and should be a prescribed smooth func- 
tion which is defined for all t > but can be singular 
when t 0. 

As we will see, only certain well-chosen functions Uq 
will be compatible with the given system of equations. 
To determine the class of compatible leading-order terms 
for a given equation, a "guess" uo must be made in a first 
step as mentioned before, often on the basis of physical 
or heuristical arguments. In a second step, certain con- 
ditions need to be checked in order to determine whether 
the chosen function uq is compatible with the equation 
in the sense that it gives rise to a well-posed singular 
initial value problem. In order to understand what this 
means, let us return to the example of the wave equation 
Eq. (TTJ) which is of second-order hyperbolic Fuchsian form 
Eq. ©, with here a = -1/2, b = 0, c = 1 and / = 0. 
The requirement that the solution u is smooth for t > 
implies that uq must be smooth in the limit t — > and 
hence the leading-order term uq is only compatible with 
the equation if it is of the form of a (truncated) Taylor 
expansion as in Eq. @. In other words, since there are 
no smooth solutions of the wave equation which become 
singular in the limit t — > 0, the leading-order term must 
be of this form. Note that if this uq consists of the first 
two terms of the Taylor expansion for example, then the 
remainder is 0{t 2 ), i.e. of higher order than uq. A partic- 
ular consequence of this is that the standard initial value 
problem of the wave equation is just a special case of a 
singular initial value problem. 

In order to determine uo beyond the case of the stan- 
dard wave equation, we will now treat the following 
"canonical set-up" which will turn out to be directly use- 
ful later in this paper. Although the following conditions 
are not always satisfied in the later discussion directly, we 
will be able to reduce our problems below to this canon- 
ical case. To this end, let us make the basic assumption 
that the right-hand side of Eq. ©, i.e. both the source- 
term function / and the second-order spatial derivative 
term, are negligible at t — in the sense that the leading- 
order behavior of the solutions to the full equation is 
driven by its principal part, only, in the now to be dis- 
cussed sense. For instance, in general relativity and the 
Einstein's vacuum field equations, the famous BKL con- 
jecture [26l429j claims that for large parts of the dynamics 
close to generic singularities, the kinetic terms in Ein- 
stein's field equations dominate the potential terms. If 
the singularity is asymptotically velocity dominated, see 
above, then this is a good description all the way to the 
singularity, and therefore our working assumption, where 
in particular all spatial derivatives are assumed to be neg- 
ligible at t = 0, is relevant. Indeed Fuchsian techniques, 
under certain analyticity assumptions, have been applied 



to asymptotically velocity dominated spacetimes before 
even beyond the Gowdy case @, 0, H3-[32| ■ In the general 
case predicted by the BKL conjecture, however, potential 
terms can lead to bounces from one kinetically dominated 
epoch to another and hence are not always negligible. As 
far as we know, Fuchsian techniques have not been ap- 
plied to such "Mixmastcr-like" singularities yet, and it is 
not clear whether this is possible. 

In any case, let us now go as far as choosing the leading- 
order term uq for the singular initial value problem to be 
a solution of the system of ordinary differential equations 
(in which now x plays the role of a parameter) obtained 
from Eq. §5§ by setting the right-hand side to zero. Since 
the coefficients of this ODE's do not depend on t, we can 
find explicit solutions easily as follows: 

fu*(x)t-< x hogt + u*4x)t-< x \ a 2 = b, 
U °~\u4x)t- x ^+u**(x)t- x ^\ a 2 ^b. [) 

Here, the smooth asymptotic data it* and be 

freely specified, and we have set 

Ai := a + \J a? — b, X2 := a — \/ a 2 — b. 

Note for later purposes that, for such a canonical two- 
term expansion, one has 

L[u ] = -t 2 c 2 dlu . 

By convention, we impose that SftAi > 5RA2, where 5R 
denotes the real part of a complex number. We refer 
to this leading-order term u$ as the canonical two- 
term expansion and the underlying argument used to 
derive it as the Fuchsian heuristics; this is motivated 
by the fact that uq is determined by Fuchsian ordinary 
differential equations. The singular initial value problem 
based on this leading-order term will be referred to as the 
standard singular initial value problem. For given 
asymptotic data, this problem is well-posed if there ex- 
ists a unique remainder w which is of ordei0 0(t a ) with 
a > — 5RA2 at t — and is smooth for t > 0. 

The question, whether the canonical two-term expan- 
sion above, or any other choice of leading-order term, 
is compatible with the given equations and gives rise 
to a well-posed singular initial value problem, depends 
strongly on properties of the right-hand side of the equa- 
tions. We cannot go into the mathematical details here 
and refer the reader to [![. The most important condi- 
tion is the following one. Consider the class of functions 
u = uq + w associated with a fixed leading-order func- 
tion uq and with all functions w that are smooth for t > 
and (together with all of their derivatives) are of order 
0(t a ) at t = for a fixed spatially-dependent function 
a. Here, it is required that a is sufficiently large so that 
w can be indeed treated as "higher order" in t. Now, 



3 See Q for the precise meaning of the symbol O in our context. 
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provided the source-term / maps each such function to 
a function which is smooth for allt > and is (including 
all derivatives) of higher order, say 0(t a+e ) for some ar- 
bitrarily small e, and provided the function L[uq] is thus 
of order 0(t a+e ), then the singular initial value problem 
can be proven to be well-posed. More precisely, this lat- 
ter condition is sufficient for well-posedness in the case of 
Fuchsian ODE's systems, only, i.e. when the right-hand 
side contains no spatial derivatives 0. In the PDE's 
set-up, there are further restrictions, in particular on the 
coefficients of the principal part, which we will, however, 
not discuss here; again see [1| for details. 

As we will discuss later, the leading-order term func- 
tion ito does not necessarily consist of only two terms as 
in Eq. (jB]). However, the free asymptotic data will in gen- 
eral consist of two free functions, at least for the class of 
second-order hyperbolic equations under consideration 
here. 

As an example, consider the Euler-Poisson- 
Darboux equation which also plays a role in general 
relativity; for example, the Gowdy equations discussed 
later (see Eq. (|TT|) . below) reduce to such an equation (in 
the so-called polarized case Q — 0, see below). For any 
real positive constant k, let 

u tt — u t = u xx . (7) 

This is a linear second-order hyperbolic Fuchsian equa- 
tion containing a single singular term. The canonical 
two-term expansion in this case reads 

u (t x) — l + u**(x)t K , k > 0, 

In [l| we showed that the singular initial value problem 
is well-posed if and only if k < 2. On the other hand 
the second-order spatial derivative term is of order 0(t 2 ) 
when k > 2, and is therefore not negligible at t = 0. Note 
that the (standard) initial value problem for the wave 
equation Eq. (TTJ is recovered in the special case K, = 1 
(and for c = 1). 



C. Approximate solutions with arbitrary accuracy 

The well-posedness theory would not be complete 
without means to actually compute the solutions under 
consideration. Our approximation scheme in [l[ precisely 
allows to compute general solutions to the singular initial 
value problem (when they exist) with arbitrary accuracy 
and can easily be implemented in practice. 

We suppose that a second-order hyperbolic Fuchsian 
system of the form Eq. ^ together with a leading-order 
term uq are given, the latter being not necessarily of the 
canonical form Eq. ((BJ) . Note that Eq. ((3]) , being a hyper- 
bolic system, gives rise to a well-posed (standard) initial 
value problem with data prescribed at any time t > 0. 



solution u 

leading-order term i 

approximation t 1 
approximation t 2 
approximation t 3 





//' 
// 







0*3*2*11 2 



t 

Figure 1. Illustration of the approximation scheme. 

This fact, together with the fact that the solution is pre- 
sumably described accurately by uo for small t > 0, can 
be exploited to approximate the solution of the singular 
initial value problem as follows; the idea is also illustrated 
in Fig. [TJ where we plot the solution u and the leading- 
order term uq schematically (as well as other functions 
to be described now). 

Observe in the figure that the leading-order function 
Mo deviates from the solution u for larger t. However, we 
can choose a small time t\ > and use the values uq and 
dtUQ at t\ as data for the (standard) initial value prob- 
lem associated with Eq. <j3j> - We solve this initial value 
problem to the future of t\ and get a function labeled 
by "approximation t\ in the plot. Since it is a solution 
of the equation which is close to the actual solution u 
at tt, it will be close to u at least on some short time 
interval. Now, the choice of t\ above was arbitrary, and 
in a next step we can pick some other time ti between 
and t\ and go through the same procedure. The result- 
ing function labeled by "approximation in the plot 
can be expected to be a better approximation of u than 
the previous function, since the deviation from u at ti is 
smaller. Eventually, we construct a whole sequence of ap- 
proximate solutions with initial times t n — > 0. If the pro- 
cedure works then this sequence approaches the actual 
solution to the singular initial value problem. The errors 
become smaller the smaller the initial time is. In fact, our 
well-posedness theory in makes use precisely of this 
approximation scheme for the existence prooQ. Therein, 
we proved convergence and also provided explicit error 
estimates. This scheme can be implemented numerically 
without difficulty, since it requires to compute a sequence 



4 Convergence and error estimates for this scheme were proven 
for linear equations, while the nonlinear situation was handled 
analytically by a further iteration, yet based on this linear case. 
Nevertheless, our numerical experiments have demonstrated that 
the scheme above does converge in the nonlinear situation. 
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of solutions to the (standard) initial value problem for 
hyperbolic equations, for which a variety of robust and 
efficient schemes are available in the literature [3^, [34[ . 

Let us make some comments about our actual numer- 
ical implementation of the above scheme; more details 
can be found in [l|, Q. First of all, recall that in order 
to compute an approximate solution in the sense above, 
we must solve Eq. (J3J) with initial data at some to > 
toward the future direction t > to- Moreover, we have to 
be able to choose to as small as needed in order to get 
an approximation as accurate as possible. The main ob- 
stacles here are the factors 1/t or 1/t 2 in the hyperbolic 
equation. We address this issue by introducing a new 
time coordinate 

t := logt. 

For instance, under this transformation, the Euler- 
Poisson-Darboux Eq. ([7]) becomes 

d 2 u — k d T u — e 2T d 2 x u = 0. 

With this coordinate choice, we have achieved that there 
is no singular term in this equation; the main price to pay, 
however, is that the singularity t = has been "shifted 
to" t — -co. Another disadvantage is that the char- 
acteristic speed of this equation (defined with respect to 
the r-coordinate) is e r and hence increases exponentially 
in time. For any explicit discretization scheme, we can 
thus expect that the so-called CFL condition^ is always 
violated after some time on. But we do not expect that 
this is a genuine problem in practice. The new time co- 
ordinate is only introduced to deal with very small times 
t. From any given larger time on, we could, if neces- 
sary, switch back to the original time coordinate t so 
that the CFL restriction above disappears. All the nu- 
merical solutions presented in this paper, however, were 
done exclusively using the time coordinate r. 

In our numerical code, we assume that a leading-order 
term uo has been fixed and, then, we write the equation 
in terms of the remainder w. The initial data, which we 
need to prescribe for each approximate solution at initial 
time To is hence 

w(tq,x)—0, d T w(To, x) = 0. 

Inspired by Kreiss et al. |3a| and by the general idea 
of the method of lines [33], we discretize the equation 
with second-order finite-differencing operators. 

D. Canonical expansions of arbitrarily high order 

One final tool remains to be presented here, that is, 
another approach to analyze the behavior of solutions 



5 The Courant-Fricdrichs-Lewy (CFL) condition for the discretiza- 
tion of hyperbolic equations with explicit schemes is discussed, 
for instance, in [33l . 



beyond the leading order. Consider the standard singular 
initial value problem associated with Eq. ^ together 
with the canonical two-term expansion Eq. ((6]). Let us 
suppose that it is well-posed, i.e. for this given uo, there 
exists a unique solution u = uo + w with remainder w 
which is smooth for t > and behaves like 0(t a ) for 
some sufficiently large a. 

Consider first the ODE's case of Eq. ([J|, i.e. 

x) + 2a(x)Dv(t, x) + b(x)v(t, x) = fo{t, x), 

where x is now regarded as a parameter and fo(t,x) is 
a given function and is smooth for t > 0. As discussed 
before, there exists a unique solution v = uo + w of this 
equation, provided suitable assumptions on fo are made 
at t = 0. At any given spatial point x, define 

(H[fo])(t, x) := r tt W f f (s, x)s a ^- 1 In - da, 
Jo s 

if a 2 (a;) 2 = b(x), and 

(H[f ])(t,x) := f f (s,x)s x ^- 1 ds 

M(x) - \ 2 {x) Jo 

f-\l(x) ft 

-77VTn / M^> Mx) -\ 
Xi{x) - \ 2 {x) Jo 

if a 2 (x) ^ b(x). Then, if these integrals converge, it can 
be checked that the solution of the singular initial value 
problem is given by 

v = u + H[fo}. (8) 

Next, consider the following ordinary Fuchsian differ- 
ential equation 

D 2 v(t, x) + 2a(x)Dv(t, x) + b(x)v(t, x) 
= f(t, x, u, tu x ,Du) + t 2 c 2 {t, x)u xx (t, x), 

which follows from Eq. (jj) for some given u = Uq + W 
where w is the remainder of u in the same sense as above. 
Here v is the new unknown. Under the conditions before, 
the singular initial value problem with leading-order term 
Mo associated with this equation for v is well-posed and 
there exists a unique solution. 

Now suppose that we start with a seed function u = uq. 
Then the singular initial value problem for Eq. ^ associ- 
ated with the leading-order term uo determines a unique 
function v, which we call u\. Then we use this func- 
tion Mi instead of u in the right-hand side of Eq. © 
and, consequently, determine a new approximate solu- 
tion v to the singular initial value problem associated 
with Eq. Q, which we call u%. Continuing inductively, 
we determine a sequence of functions (u n ), each of these 
being computed using the integral expressions above and 
having uo as their leading-order term. It was observed in 
[2, 0, [I] that Uj+i - uj is 0(t^) at t = with {j3j) some 
monotonically increasing sequence of constants. Hence 
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Uj can be interpreted as a formal expansion of the solu- 
tion at t = whose order in t increases with j. Moreover, 
it turns out that the residual, which is obtained when 
Uj is plugged into the full equation Eq. (@}, behaves like 
a positive power in t at t = which increases with j. 

We note that one can find conditions for which this 
sequence (u n ) actually converges to the solution of the 
singular initial value problem in general only in the case 
of Fuchsian ordinary differential equations. In the PDE's 
case, the sequence (u n ) does not converge in general. 
Even though it may not converge, it is still useful for 
a qualitative study of general solutions. If one is rather 
interested in quantitative statements and in actual con- 
vergence, then the approximation procedure in the pre- 
vious subsection must be used. 



Let us now proceed with a heuristic discussion of the 
equations in order to motivate the choice of the leading- 
order term for the singular initial value problem. Based 
on extensive numerical experiments [TEl Il6l |36| and the 
ideas underlying the (already mentioned) BKL conjec- 
ture, it is conjectured that the spatial derivatives of solu- 
tions (P,Q) to (fTTj) becomes negligible as one approaches 
the singularity and, hence, (P, Q) approach a solution of 
the ordinary differential equations 



Ptt + — = e Q t , Qtt 



Qt 
t 



= -2R 



These equations are referred to as the velocity term 
dominated (VTD) equations. They admit solutions 
that are given explicitly by 



III. GEODESICS IN GOWDY SPACETIMES 

A. Background material 

Before we apply the general theory, we provide some 
background material about Gowdy spacetimes. These 
are spacetimes with two commuting spatial Killing vec- 
tor fields, for which one can thus introduce coordinates 
(*, x, y, z) so that the functions y, z are aligned with the 
symmetries. The so-called Gowdy metric can then be 
written in the form 



g = J_ e V 2 (-dt 2 +dX 2 ) 

+ t(e p (dy + Qdz) 2 4 



(10) 



with t > 0. It therefore depends on three coefficients 
P = P{t.x), Q = Q{t,x), and A = A(t,x). We as- 
sume that these functions are periodic with respect to x. 
Clearly, the metric is singular (in some sense) at * = 0. 
We recall that a Gowdy metric is said to be polarized 
if the function Q vanishes. 

Einstein's vacuum equations imply the following 
second-order wave equations for P, Q: 



Pu + j-P x , 
Qtt H — - — Qxx 



e 2F (Q 2 t-Ql), 
—2(PtQt — PxQa 



(11) 



which are decoupled from the wave equation satisfied by 
the third coefficient A: 



A, 



A, 



Pi - Pt + e^{Qi - Qt 



(12) 



Moreover, the Einstein equations imply the following 
constraint equations: 



At 



2* (P x P t + e 2P Q x Q 



t (Pi 



te 



Pi 



Q 



(13) 



Eqs. (fTTj) represent the essential set of Einstein's field 
equations for Gowdy spacetimes; cf. |16| for further de- 
tails. We refer to them as the Gowdy equations. 



P(t,a:) = log {at k (l + C 2 t- 2k )), 
(t~ 2k 

Q(tlX) a(l + C 2 i~ 2fc )' 



(14) 



where x plays simply the role of a parameter and a > 0, 
C, £, k are arbitrary 27r-periodic functions of x. We as- 
sume in the following that k > and ( ^ 0. Based on 
the above formulas, it is a simple matter to determine 
the expansion of the function P near t = 0, that is, 

Pit x) 

lim v ' ' = limtP t (t,x) = -k, 

t^Q log* t-X) ^ ' 

and 

lim (P(t, x) + k(x) log*) = log(aC 2 ). 
Similarly, for the function Q we obtain 
fmQ(t,x)=£-± 
lim*- 2 ' 5 (Q(* >a; )-£ + i-) =-L. 



We thus arrive at the expansions 

P(t, x) =-k(x) log * + P** (x) + . . . , 
Q(t,x) = Q*(x)+t 2k WQ**(x) + ..., 



(15) 



in which k, P**, Q*, Q** are functions of x. In general, 
P blows up to +00 when one approaches the singularity, 
while Q remains bounded. Note that if the assumption 
k > is dropped in Eq. (TT4")) . then it turns out that we 
must substitute k by \k\ in most of the previous expres- 
sions. In other words, without loss of generality we can 
assume k > for the leading-order expansion Eq. (|15l) : 
we will ignore the exceptional case k = in most of the 
following discussion. 

It was shown in [l] that this leading-order term 
Eq. (|15[) is the canonical two-term expansion of the 
Gowdy equations. Indeed, we found that the singular ini- 
tial value problem is well-posed as long as k is a function 
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with values in the interval (0, 1). If d x Q*{x) — at points 
where k > 1, then the singular initial value problem is 
also well-posed. These restrictions are related to the for- 
mation of spikes, investigated earlier in [IH, [H, [36[ . For 
asymptotic data Q» = Q*„ = 0, the corresponding solu- 
tion is polarized and then the function k(x) may take all 
values in R. Given a solution of Eq. (ITT1) with a leading- 
order term Eq. (|15[) satisfying these restrictions, then also 
Eq. (1121) can be solved as a singular initial value problem 
with leading-order term 



A(t,x) = A*(x)\ogt + A**(x) + 

where 

A, (a;) =fc 2 (x), 

A„(x)=A + 2 f (-Pi,+2e 2P "Q^Q',)kdx. 
Jo 



(16) 



(17) 



Here, a prime denotes the derivative with respect to x. 
Note that periodicity in space hence implies the following 
further restriction on the asymptotic data 



/ {-PL 

Jo 



2e 2F **Q**Q'J kdx = 0, 



and the constraints Eqs. (|13[) are then satisfied for all 
t > 0. It was demonstrated in (3(| that solutions of the 
Gowdy equations which have the leading-order behavior 
as above approach certain Kasner solutions at t — 0, with 
in general different parameters along different timelines 
toward the singularity. 



B. Timelike and null geodesies in Gowdy 
spacetimes 

Now we know how to construct Gowdy solutions 
of Einstein's vacuum field equations with a prescribed 
asymptotic behavior, and we can assume that such a 
spacetime is given. We henceforth study freely falling ob- 
servers and null geodesies in a vicinity of the time t = 0. 

A geodesic curve 7 A '(s) is a solution of 



~ds T 



r A(7(*))^(s)^(s) 



0. 



where s is the affine parameter which, in the time- 
like case, corresponds to the proper time of the ob- 
server traveling along the geodesic. The objects „ 
are the Christoffel symbols of the Gowdy spacetimes. 
For the later discussion it will be more convenient to 
parametrize the geodesies with respect to the coordinate 
time t(s) = 7 (s). Hence from now on we will consider 
the functions 7 P as functions of t and a dot will denote the 
derivative with respect to t. We will assume future point- 
ing causal geodesies with dt(s)/ds > 0. The geodesic 
equation becomes 



and our main aim in the following is to study the singu- 
lar initial value problem for this equation. Since it is a 
system of ODE's, the analysis in the part of this paper 
simplifies. 

Note that if 7 M (t) is a solution of Eq. (jTSJ) we can return 
to using the affine parameter s(t), determined from the 
equation 



const 



£) (ffoo+^fV), (19) 



where a, /3 = 1,2, 3 are the spatial coordinate indices of 
the metric g^ v of the form Eq. (flUl) . The constant in this 
equation is positive for a timelike geodesic and zero for a 
null geodesic, and may be specified freely. 



C. Orthogonal causal geodesies 

Our strategy will be to increase the level of complexity 
of our problems systematically step by step. We start by 
discussing causal geodesies which are orthogonal to the 
Gowdy symmetry orbits, i.e. 7 y = r y z = 0. We refer to 
those as orthogonal geodesies. We find easily that the 



conditions 7* 



7 Z = at some initial time t > implies 
for all t > 0. We write x(t) instead of "/ x (t) 
for the only remaining non-trivial component, and the 
geodesic equation (TT5|) reduces to 



- ^d x A(t,x{t))x 2 (t) + ^d x A(t,x(t)) = 0. 



(20) 



We now study Eq. (|2U|) as a singular initial value prob- 
lem. The first step is to "guess" the leading-order be- 
havior of the geodesies at t = 0. Recall that the Gowdy 
solutions considered above have the property that spa- 
tial inhomogeneities close to t — are insignificant. This 
suggests that geodesies should behave, to leading order, 
like in spatially homogeneous (Bianchi I) spacetimes. We 
will find in the course of the discussion that this does lead 
to the correct leading-order term for the geodesies in the 
general Gowdy case. 

a. Step 1. Geodesies in Kasner spacetimes. Kas- 
ner spacetimes are particular spatially homogeneous so- 
lutions of the Einstein's vacuum field equations [37| . 
The Kasner metric can be brought to the Gowdy form 
Eq. (flU)) by requiring that 



P(t,x) = -klogt + P**, Q(t,x)=0, 
A(t,x) = k 2 logt + A^, 



(21) 



r + (r„% - r v ° P r) rr = o, 



(18) 



where fc, P** and A** are arbitrary real constants. Hence 
Kasner solutions can be considered as solutions to the 
singular initial value problem for the Gowdy equations 
above with asymptotic data constants fc, A**, cf. 

Eq. (TTT]) . Note that k is allowed to be any value in R here. 
Moreover, the constants P** and A** can be transformed 
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to zero by suitable coordinate transformations. We thus 
have 

g = t 1 ^ (-dt 2 + dx 2 ) + t x - k dy 2 + t 1+k dz 2 

and, in the notation made in [37j for the Kasner expo- 
nents, 

Pl = (k 2 -l)/(k 2 + 3), 
P2 = 2(1 -k)/(k 2 + 3), 
P3 = 2(1 + k)/(k 2 + 3), 

so that the three flat Kasner cases are realized by k = 1, 
k = — 1, and |fc| — > oo, respectively. 

For any Kasner spacetime, Eq. (I2TJ)) reduces to 



1-k 2 

At 



(x(t)-x 3 (t)) =0. 



This equation can be integrated once 

XlT](t) 



±(t) = F(r,(t)) 



v / T+WW' 



(22) 



(23) 



where 



n(t) = t* 



and X\ is an integration constant. Since F is smooth in 
r\ for all r\ € R, we can compute its Taylor expansion at 
7/ = 0. If \k\ < 1, this yields the leading-order behavior 
of x(t) at t = 0. If |fc| > 1, we set e := i]" 1 and consider 
the function 



G(e) :=F(l/e)= Xl 



1 



for e > 0. This function can be extended as a smooth 
function to e < in only one way and then we are allowed 
to compute its Taylor expansion at e = 0. If |fc| = 1, the 
exact solution of Eq. (|22|) is x(t) = const. Then, after a 
second integration, the leading terms of the expansions 
for x(t) at t = can be computed: 



x, + ^ Xl t^- k2 ) + 

Xq 



x{t) = < 



x + sign(a;i) t 



sign(xi) _ t i {1+ k 2 ) 



x\(l 



\k\ < 1, 
\k\ = 1, 

\k\ > 1, 



(24) 



where xq is another integration constant and x\ ^ 0. We 
have introduced the sign function 



sign(a;) 



Observe that the case x\ — corresponds to having 
x(t) = xq for all three cases. 




The dots represent terms with t-powers n(\ — k 2 )/A+l 
for n = 2, 3, . . . in the first case, and n(k 2 — l)/2 + 1 in 
the third case. The expression for \k\ = 1 is exact. Given 
a sequence of Kasner solutions all with, say, |fc| < 1 (and 
similarly when |fc| > 1) for which \k\ approaches 1, the ex- 
ponents in the i-powers for all of the terms in the infinite 
series above approach to 1. Th e scries given by the sum of 
these terms converges to x\j \fl + x\. This demonstrates 
the important fact that the closer we choose |fc| to the 
value 1 , the larger all higher-order terms in the series be- 
come which, consequently, become no longer negligible. 

For all values of xq^xi G R the orthogonal geodesies 
above are timelikc. Null rays, which are given by x(t) = 
±1, are obtained formally in the limit x\ — > ±oo. Now, 
Fig. [U shows the leading-order behavior of orthogonal 
geodesies close to t = for the three cases \k\ < 1 and 
|fc| > 1. In the first case, the geodesies approach a curve 
at t = which is at rest with respect to the coordinate 
system, while in the third case, the geodesies approach 
a null geodesic. Recall that the case |A;| = 1 corresponds 
to flat Kasner solutions and hence the curves are straight 
lines. 

b. Step 2. Fuchsian analysis for the Kasner case. 
Instead of deriving Eq. (|2"4")l as above, let us now ana- 
lyze Eq. (|22p by means of the Fuchsian heuristics. We 
multiply Eq. (|22[) with t 2 and write the time derivatives 
by means of the operator D = td t . 



D 2 x--(5-k 2 )Dx 



1 



At 2 



(Dx) 



(25) 



This equation is of second-order Fuchsian type. For the 
standard singular initial value problem for this equation, 
where we interpret the left-hand side as the principal 
part and the right-hand side as the source-term, we seek 
solutions with leading-order behavior 



x(t) 



+ x^ti {5 k2) + 



(26) 



at t = for arbitrary prescribed constants a;*,x** € 
R if k 2 < 5. For arbitrary remainders w — 0(t a ) 
with a > j(5 — k 2 ), the source-term is t~ 2 (Dx) 3 = 

0(ti( 5 ~ k2 ^~ 2 ) in general or identically zero if |fc| = 1. 
Hence, in general, the source-term is negligible if the fol- 
lowing is positive 

(3(5 - k 2 )/A - 2) - (5 - k 2 )/A = (1 - fc 2 )/2, 

since then, the quantity a can be chosen between (5 — 
k 2 )/A and (3(5 - k 2 ) /A - 2) for all derivatives. Therefore 
the singular initial value problem for Eq. (|2"5|) associated 
with Eq. (121)1) is well-posed if |fc| < 1. It is certainly 
also well-posed for |fc| = 1 since then the source-term 
is identically zero. Indeed, for |fc| < 1, the solutions of 
the singular initial value problem agree, for appropriate 
choices of a;* and with the previous explicit solutions 
with expansions Eq. 

For |fc| > 1, the singular initial value problem for 
Eq. (1231) associated with Eq. (|2"6")l is not well-posed, since 
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Figure 2. Orthogonal geodesies in the homogeneous polarized 
case for k = 0.5, k = 1.8. 



the source-term is not negligible at t = 0. Note that 
Eq. ([2"B)) is also not compatible with Eq. (pM)) in this case. 
How can we proceed? The expansion Eq. does not 
have the form of a canonical two-term expansion, due to 
the additional term ±t. Suppose that the right-hand side 
of Eq. (|2"5)) is not negligible but rather, say, of the same 
order in t at t — as the left side. We make the ansatz 
x = x* + t^c + . . . where x„, (3 and c are unknown con- 
stants. Then the equation implies that f3 = 1 and that 
c — 1, independently of k. This explains the presence 
of the term ±t in Eq. (|24)) . It suggests that we should 
introduce new variables 



f±(t) :=x(t)Tt, 



(27) 



for the cases x\ > and x\ < in Eq. (|2~4")) . In terms of 
those, Eq. (|25l reads 



(28) 



D 2 f ± -1(1 + e)Df ± = T [Df ± f 



1 



At 2 



{Df±? 



This is a second-order Fuchsian equation whose canonical 
two-term expansion is 



(29) 



at t = with arbitrary asymptotic data G M. 

This clearly looks promising, cf. Eq. ((24)) . We have to 
check whether the source-term, i.e. the right-hand side 
of this equation which consists of two terms is negligible 
for arbitrary remainders w — 0(t a ) with a > \{\ + k 2 ). 
For the first term, we get 



2 = D(f( 1+k )"! 



), 



and hence the following must be positive 

(l + fc 2 )-l-i(l + fc 2 ) = i(fc 2 -l). 

This is true for |fc| > 1. We find that also the second 
term is negligible for \k\ > 1. Hence the singular initial 
value problem for Eq. (j28|) associated with Eq. (l29l) is 
well-posed for \k\ > 1. The case = corresponds 
to null geodesies. We get timelike geodesies if we choose 
/ + for x** < 0, and, if we choose /_ for x* r > 0. The 
case x\ — before now formally corresponds to the limit 
— > oo. 

The leading-order terms for the Fuchsian analysis are 
therefore consistent with Eq. ((24)) in the Kasner case. A 
particular result of our efforts so far is that all causal or- 
thogonal geodesies in Kasner spacetimes have limit points 
at t = 0, and hence "start" from points on the singularity 
which can be prescribed freely. We will see that this is 
not always the case for non-orthogonal geodesies. 

c. Step 3. The general equation for orthogonal 
geodesies. Now we proceed with the singular initial 
value problem for the general Gowdy inhomogeneous 
(unpolarized) case for orthogonal causal geodesies. The 
leading-order behavior which we have identified in the 
Kasner case before will be taken as the guess for the 
leading-order behavior here. As before, we consider the 
function A as a given solution of Eq. ([T2")) compatible with 
Eq. (TTTj)) and (fTT)) at t = where, in particular, k(x) is 
now any given function with the restrictions discussed 
earlier; we assume here that it is positive. Again, we 
rewrite Eq. (|20[) as a second-order Fuchsian equation 

1 



D 2 x(t) - -(5 - k\x*))Dx(t) 

= ^Dx(t)(k 2 {x*) - DA(t,x(t))) 



1 - DA(t,x(t)) 
it 2 



(30) 



[Dx{t)f 
t 2 )d x A(t,x(t)). 



Here, DA(t,x(t)) means the partial derivative of A with 
respect to its first argument evaluated at (t, x(t)) and 
multiplied by t. Similarly, d x A(t, x(t)) means the par- 
tial derivative with respect to the second argument eval- 
uated at (t,x(t)). Note that we have added the term 
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k 2 {x*)Dx(t)/A to both sides of Eq. the significance 
of this term for the Fuchsian analysis becomes clear in a 
moment. 

Let us suppose first that k(x#) < 1 for some x* £ 
R. The additional term k 2 (x*)Dx(t) / A in both sides of 
the equation has the consequence that the canonical two- 
term expansion coincides with that of the Kasner case 
Eq. ([24|. Hence we shall consider the singular initial 
value problem associated with the leading-order term 



xi^t^j — X% ~\~~ X-^*t ^ ^ ^ ^ — I— 



(31) 



for arbitrary x*,x** € K. It follows that the source-term 
is negligible. In order to show this, we note that 

k 2 (x,)-DA(t,x(t)) = 0(t a ), 

for some a > 0, and 

d x A(t,x(t)) = 0(logt). 

Hence, as expected, the singular initial value problem for 
Eq. ([50)1 associated with Eq. (T5TT) is well-posed if k(x*) < 
1. The same conclusion follows directly for fc(x») = 1. 
In order to study the case A; (a;*) > 1, we make the same 
ansatz as in Eq. (|27|) . The resulting equation can be 
analyzed in exactly the same way in well-posedness can 
be shown. We omit the details for this case. 

So, causal future directed orthogonal geodesies have 
the asymptotic behavior given by Eq. with k — 

k(x*). Hence we have shown that spatial inhomogeneities 
as well as the polarization of the given Gowdy spacetime 
do not play a significant role for the leading-order dynam- 
ics of observers close to the singular time t = 0. Note that 
Gowdy solutions where k takes values only in the interval 
(0, 1) have a curvature singularity at t = 0. In particular, 
this singularity is asymptotically velocity dominated as 
mentioned before and therefore consistent with the BKL 
conjecture. 



D. Non-orthogonal causal geodesies 

Next we consider geodesic curves 7 A '(i) which are not 
necessarily orthogonal to the Gowdy group orbits. The 
full discussion of this problem would, however, be too 
lengthy for this paper since all three spatial coordinate 
components of 7^ would then be non-trivial in general. 
This is why we restrict to the polarized Gowdy case 
Q = 0. In this case, the functions 7 Z and 7 y decouple in 
Eq. (|T5| . We assume that j v = and 7 Z is non-trivial. 
There is no loss of generality since one can always switch 
to the case 7 Z = and a non-trivial 7 y by using the 
transformation k h-> — k. The latter mapping just inter- 
changes the two Killing vector fields and, hence, is an 
isometry of the spacetime. We will be particularly inter- 
ested in the case < \k\ < 1, but also in the borderline 
case \k\ = 1, due to its relevance for the formation of 
Cauchy horizons, as will be explained below. 



As before, we write x(t) and z{t) for the now two non- 
trivial coordinate components of the geodesic curves. We 
get the following coupled system of second-order Fuch- 
sian equations 

D 2 x(t) - -(5 - k 2 (x*))Dx(t) 

= -t 2 F 1 {t 1 x{t)) 

+ F 2 (t, x{t))Dx(t) + Fi(t, x{t))(Dx{t)) 2 

+ F 3 (t, x{t)){Dz{t)) 2 + F A (t, x{t)){Dx{t)f 

+ F 5 (t,x(t))(Dz(t)) 2 Dx(t), (32) 

D 2 z(t) - ^{k 2 (x r ) - 4fc(x*) - l)Dz{t) 

= Gi(<, x(t))Dz(t) + G 2 (t, x{t))Dz{t)Dx{t) 
+ F b {t,x(t))(Dz(t)f 
+ F 4 (t,a;(t))L>z(t)(L>a;(t)) 2 , 



with 



Fi(t,x) = -8 x A(t,x), 

F 2 {t,x) = - A {k 2 { Xar )- DA(t,x)), 

F 3 (t,x) = -i e -( p (^)+ A (*-)/ 2 )t 3 / 2 9 :c P(i,x), 



Fi(t,x) = - 



1 - DA(t,x) 
it 5 ■ 



F 5 (t, X) = I e -(*<t.*)+A(t,*}/2) t -l/3 (1 _ Dp(t) 



and 



Gx(t, x) = + DP(t, x)-~ (k 2 (x*) - DA(t, x)) , 

G 2 (t,x) =d x P{t,x) + ^d x A(t,x). 

Note that we have added certain terms to both sides of 
the equation as before; the second term on the right- 
hand side of the equation of x would not satisfy the de- 
cay condition for well-posedness without the additional 
term k 2 (x*)Dx. The same is true for the additional term 
— j(k 2 (x*) — 4:k(x*))Dz. However, it turns out that the 
standard singular initial value problem for this system 



z(i) = z, + z^(fe 2 (^)-4^-)-i) + 



(33) 



is still not well-posed. Consider for instance the term 
F${t, x{t)){Dz(t)Y which does not have the correct be- 
havior at t = 0. This suggests that Eq. (|3"3"|) is not the 
correct leading-order behavior. Let us hence step back 
and think about the dominant physical effects in the same 
way as we did for the case of orthogonal geodesies. 

d. Step 1. A decoupled equation for \k\ < 1. We 
have realized before that inhomogeneities do not con- 
tribute to the leading-order term of orthogonal geodesies. 
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This suggests that the Kasner case given by Eq. (f^Tj) with 
constant parameters k, P** and A** is also relevant here. 
However, it turns out that the standard singular initial 
value problem remains ill-posed when =/= 1. Moreover, 
the system is still too complicated to find explicit solu- 
tions and hence we cannot proceed as for instance for 
Eq. ([231). 

So let us simplify the equations even more. We can try 
to neglect the coupling between both equations Eq. (1521 . 
For this we consider the equation for z(t) with x(t) 
( i (of course the equation for x with vanishing z is the 
geodesic equation in the orthogonal case already treated) : 



D 2 z{t) 



1 



(fc 2 -4k- l)Dz(t) 



= I e -(^.+A«/2);-(l-fc) 2 /2 (1 + k){Dz[t)) \ 



(34) 



The standard singular initial value problem is still ill- 
posed. However, we are able to find explicit solutions for 
z(t), namely 



z(t) = ±F(?7(t))t-i (3 - fc) ( fc+1 ) or z(t) = 0, 



with 



and 



F(rj) 



e i 



(P«+A,„/2) 



i(P«,+A»»/2) 



2_Zi7? 



n(t) = t 1+k , 

for arbitrary z\ € K (provided 77 is sufficiently small) . As 
before, it turns out that F is a smooth function in 77 if 
77 is small, and its Taylor expansion at 77 — henceforth 
yields an expansion at t — 0, as long as fc > — 1. After 
an integration, we then find that 



z(t) =zq ± e 



±Zi 



i(P«+A«/2). 



5 + 2k + k 2 



t~ 



(1-fc) 2 

i(5+2fc+fc 



or 



Z(t) = Zq, 

both with another free parameter zo £ i. We can see 
clearly the reason why the standard singular initial value 
problem Eq. ([33)) has failed. Namely, this expression is 
not of the form of a canonical two-term expansion be- 
cause the second term is of lower order at t = than the 
third term if k > — 1 . Similar to the case of orthogonal 
geodesies, cf. Eqs. ([27)) and (|28|) . we define 



f±(t) :=z(t) T e^ P ~+ A ~^- 



rtH 1 -^. (35) 



(1-fc) 2 

Indeed, this ansatz turns Eq. (|34p into a second-order 
Fuchsian equation for f± with a canonical two-term ex- 
pansion 

f ± (t)=z*+z^ 5 + 2k+k2 ) + ... 



We find easily that the corresponding singular initial 
value problem is well-posed for \k\ < 1. 

e. Step 2. Kasner case for \k\ < 1. We have consid- 
ered Eq. ([Ml) , in order to get an insight about the behav- 
ior of geodesies. Now we check that the same leading- 
order term for z(t) as given by Eq. (|3"5|) with f±(t) de- 
fined in Eq. ([35]) is consistent with the fully coupled set 
of equation Eqs. (j52"j) for the Kasner case with \k\ < 1. It 
is important to note that the ansatz Eq. (j55"j) for z does 
also change the principal part of the equation for x due to 
terms which have originally been part of the source-term. 
The equations become 

D 2 x(t) - -(7 + 2k - k 2 )Dx(t) = 

D 2 f±(t) - ^(5 + 2k + k 2 )Df ± (t) =..., 

where we do not write the lengthy source-terms. The 
standard singular initial value problem associated with 
the leading-order terms 

f ± (t) = z* + z^ 5+2k+k ^ + ..., 

turns out to be well-posed if \k\ < 1, i.e. all terms in the 
source-term are negligible. It is easy to see from Eq. ([T9]) 
that the geodesies here are timelike if we choose /+ for 
z** < 0, and if we choose /_ for z** > 0. The other 
solution given by z** = is ignored in the following, since 
it corresponds to the orthogonal case with z(t) = const 
which has been investigated before. 

It is very interesting to note that the leading-order 
dynamics in the x-direction is significantly different than 
in the orthogonal case, cf. Eq. (j24|) . and there does not 
seem to be a simple way to describe the transition from 
the non-orthogonal to the orthogonal case at the level of 
the leading-order terms. All the geodesies above have a 
limit at t = in the same way as in the orthogonal case. 
We can parametrize the leading-order terms of timelike 
geodesies as follows: 



x(t)=x,+x^ +2k - k ^ + ..., 
z(t) =z* - sign(z«)e' (P " +A " /2) 
+ z„t^ +2k+k '^ + . . . 



(1-fc): 



/. Step 3. Kasner case for \k\ = 1. In the flat Kas- 
ner cases |fc| = 1, we can solve the fully coupled system 
Eq. ([32)) explicitly for ( x' (t), z' (t)). From this, one gets 
the following leading-order behavior 



x(t) = 



z{t) 



xq + x\t, k = — 1, 

X + \xxt 2 + . . . , fc = +1, 

Zq + Zit, k = 

Zo±e i(P,.+AW2) los ; 



1 9 
± -Zit 2 

2 



fc 



(36) 
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for arbitrary parameters Xo, %i, zq, z\ G R. The expres- 
sions for the case k = — 1 are exact (hence no dots), 
while the dots in the case k — 1 represent terms with 
even powers in t. In the case k = 1, the leading-order 
behavior is not of the form of a canonical two-term ex- 
pansion because of the logarithmic term. In the same 
way as before, we define a new quantity for z. Then, it 
is straightforward to show that the singular initial value 
problem associated with the leading-order term given by 
Eq. ([31))) is well-posed for \k\ = 1. 

The case k = 1 is outstanding so far because the 
geodesies have no limit point at t = 0. Rather, the 
geodesies "swirl" toward the t = 0-surface in the z- 
direction. We can also conclude that this "swirl behav- 
ior" also occurs for k = — 1 for geodesies with non-trivial 
dynamics in the y-direction. 

g. Step 4- General inhomogeneous polarized case. 
The final step is to insure that the leading-order behav- 
ior, which we have identified in the homogeneous case 
now, is also valid in the inhomogeneous polarized case, 
at any point (a;*, z*) at t = 0, with k substituted by 
k(x„). For most of the terms in the now fully general 
source-term of Eqs. (f3"2")) . negligibility follows easily. How- 
ever, for some of them we need to know the leading-order 
behavior of the remainders w\ and W3 of P and A, re- 
spectively. From the theory of higher-order canonical 
expansions discussed before, it follows 

w^x) = \t 2 (k"(x)(l - log*) + PiUx)) + 0(t 4 )- 

We are allowed to take the derivative of this expansion 
and hence obtain that Dw\ = 0(t 2 \ogt). Now, Eqs. (fT3"|). 
together with Eq. (ITrJ)) and (fTT|) . allows to derive that 
W3,Dw3 — 0(t 2 log 2 t). This is sufficient to conclude that 
the singular initial value problem with the above leading- 
order term is well-posed in the inhomogeneous case if 
|fc(a;*)| < 1. Finally, in the two cases k(x*) = ±1, we 
find that the source-term is negligible and hence that the 
singular initial value problem associated with Eq. f|36[) is 
well-posed if and only if 

PU(x*) = fc'(x*) = 0, Pi'M = k"(x*) = 0. 

In summary, we have found again that, as long as k{x*) < 
1, the leading-order behavior of the geodesies is not ef- 
fected by inhomogeneities. 

IV. NUMERICAL EXPERIMENTS 

We will investigate the following situation now. As 
we have already mentioned, the (standard) initial value 
problem for Einstein's vacuum field equations is well- 
posed. It is also known from the classical work [38| that 
every vacuum initial data set generates a unique, so- 
called maximal globally hyperbolic development. 
This is the maximally extended spacetime which is, 
roughly speaking, uniquely determined by the initial data 



and within which causality holds. The example of the 
Taub solution [3^, H(| shows, however, that there are ex- 
amples of spacetimes where this maximal development 
can be extended in various ways so that the correspond- 
ing initial data do not determine the fate of all observers 
in the spacetime uniquely. Moreover, in such a space- 
time, there exist closed causal curves and causality is 
thus violated. These undesired properties have caused 
an ongoing debate in the literature whether such patho- 
logical phenomena are typical features of Einstein's the- 
ory of gravity - in which case it could not be considered 
as a "complete" physical theory - or whether such phe- 
nomena only occur under very special and non-generic 
conditions — for instance, the high symmetry of the Taub 
solutions was pointed out. 

In this context, an interesting and intensively debated 
hypothesis is provided by the strong cosmic censor- 
ship conjecturqj whose widely accepted formulation 
was proposed in [4l|, based on ideas in [42] going back 
to Penrose's pioneering work [43| . This conjecture states 
that spacetimes such as the Taub spacetimes are non- 
generic and can occur only for very special sets of initial 
data. We note that this conjecture has been proven so 
far only in special situations (covering spacetimes with 
symmetries, only). 

Now, if the maximal globally hyperbolic development 
of a given initial data set can be extended, as it is the case 
for the Taub solutions, then its (smooth, say) boundary, 
which is a so-called Cauchy horizon, signals a break- 
down of determinism for Einstein's theory, at least for 
the particular choice of initial data under consideration. 

In order to get some insights about the possible patho- 
logical behaviors arising in Einstein's theory, we propose 
here to precisely investigate solutions with Cauchy hori- 
zons; see also [44| for earlier investigations. In the Gowdy 
case under consideration here, such spacetimes can be 
generated with the singular initial value problem, dis- 
cussed in the first part of this paper, provided a special 
choice of asymptotic data functions fc, P**, Q* and Q»» 
is made. 

It turns out that a solution of this singular initial value 
problem has a curvature singularity at t = at all those 
values x where k =/= 1. On the other hand, when k = 1 on 
an open interval, then there is instead a Cauchy horizon 
there at t = and curvature is finite. We note that 
the singular initial value problem is possibly the only 
reliable way to construct a large class of Gowdy solutions 
with Cauchy horizons. For the (standard) initial value 
problem, it is in general not known which choice of initial 
data at some t > corresponds to a Cauchy horizon at 
t = on the one hand. On the other hand, even if we 
knew the data, the numerical evolution would in general 
be highly instable. 



The strong cosmic censorship conjecture should not be confused 
with the so-called weak cosmic censorship conjecture about 
the existence of event horizons in the context of black holes. 
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x = 2tt 



Figure 3. Schematic diagram for the solution considered here. 



Our particular class of asymptotic data of interest in 
the following is given as follows: 



k{x) = 



1- 



x £ [it, 2tt], 



l-e-V-e-i/^-), xe(0,7r), 
P**(x) = 1/2, Q*(x)=Q, 



0. 



-l/x p -l/(ir-») 



a; £ [7T, 27r] 
a; G (0, 7r), 



A*(x) = k 2 (x), A**(x) 



This is the same class of data considered earlier in [l|, Q ■ 
Observe that the asymptotic data and, hence, the corre- 
sponding solution for all t > are smooth but not ana- 
lytic. The solution has a curvature singularity at t = 
for x £ (0, 7r) and a smooth Cauchy horizon at t = 
for a; € (71", 27r). The spacetime looks schematically as 
in Fig. [21 This figure shows two relevant regions with 
respect to the t- and x-coordinates. Recall that the re- 
maining two coordinates y and z correspond to symmetry 
directions so that each point in this figure actually rep- 
resents a spacelike 2-torus perpendicular to the pictured 
plane. Also recall that the rr-direction is 27r-periodic. 

For the following discussion it will be convenient to 
reverse the time orientation, i.e. to think of the surface 
t = as the future boundary of the spacetime t > 0. Due 
to the particular structure of the Gowdy metric Eq. ([TU| , 
the light rays in the t, x-direction travel along 45° straight 
lines. The left-hand dark region in Fig. [3] is hence that 
part of the spacetime from which no observer nor light 
ray can escape; all future directed causal curves in there 
will hit the gravitational singularity. Hence this region 
shares many properties with a black hole in the asymp- 
totically flat case, and we will call this the "black hole- like 
region" . On the other hand, the light shaded region on 
the right is that part of the spacetime in which all ob- 
servers and light rays will eventually cross the Cauchy 
horizon. Since the Cauchy horizon is a perfectly smooth 
surface and the gravitational field is finite there, all ob- 
servers will be able to continue beyond it. Indeed, for 
the special choice of data above the whole light shaded 
region is (locally) a part of Minkowski space, i.e. there 
is no gravitational field at all. Note, however, that the 
Cauchy horizon is generated by closed null curves and 
hence causality breaks down. Nevertheless, we call this 




Figure 4. Contour plot of the absolute value of the 
Kretschmann scalar. 



region the "Minkowski region". So, the singular initial 
value problem allows us to construct solutions with quite 
peculiar properties. We expect that the part of the space- 
time limited by the black hole-like and Minkowski regions 
to be of particular interest. 

We discuss now our numerical results of this space- 
time. We have demonstrated convergence of our numeri- 
cal scheme in [l| already. In the following, we present one 
approximate solution which was computed with initial 
time at t = —10, i.e. t ~ 5 • 10 -5 , with spatial resolution 
Ax = 1CT 3 and time resolution At = 0.0025. In Fig. H 
we show a contour plot of the Kretschmann scalar K of 
the spacetime, that is, the squared norm of the Riemann 
curvature: 

This is a non-trivial curvature invariant which can be 
considered as a measure of the gravitational strength. As 
expected, K becomes very large close to t — inside the 
black hole-like region; indeed the numerical values reach 
10 15 very close to t = 0. The curvature is zero insight 
the Minkowski region with a maximal numerical error 
10~ 2 close to the upper tip. We consider this accuracy 
as sufficient for the discussion here; it would, however, be 
straightforward to increase the numerical accuracy even 
further by choosing a higher resolution or a smaller initial 
time of the approximate solution. Note that the fact that 
the curves in Fig. 0] and the following plots are slightly 
fuzzy is not caused by a lack of numerical accuracy of 
our code, but rather by the algorithm used by gnuplot 
to find contour lines. The plot reveals the extremely 
interesting dynamics here. For example note that there 
is an almost, but apparently not exactly flat region in the 
upper left corner of the figure, while there is a region in 
the upper right corner where the gravitational field seems 
to be concentrated before it collapses to the future inside 
the black hole-like region. 

Motivated by our schematic understanding of the 
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Figure 5. Causal geodesies from the intermediate point. 

spacetime, we consider the following simple "experiment" 
now. Since the gravitational field collapses inside the 
black hole-like region and hence becomes increasingly 
strong there on the one hand, and since the field van- 
ishes in the Minkowski region on the other hand, we ex- 
pect that freely falling observers traveling towards t = 
from t > are attracted to the left-hand x-direction. 
It can be expected that these observers can only evade 
the black hole-like region if they are boosted sufficiently 
strongly to the right (but not too strongly due to peri- 
odicity in x). The borderline case when the observers hit 
the intermediate point, given by t — and x — tt, and 
therefore "just miss" the singularity is hence particularly 
interesting. 

Our singular initial value problem for orthogonal 
geodesies above allows us to compute all causal orthog- 
onal geodesies that end in this point. The numerical re- 
sult is shown in Fig. [5l There we plot several orthogonal 
geodesies of this kind and we zoom into the region around 
the intermediate point. (The curve labeled by "a=0.000" 
will be explained shortly.) These numerical results have 
been obtained by solving the singular initial value prob- 
lem for the orthogonal geodesic equation numerically on 
the numerically generated background spacetime. For 
this we use second-order interpolation in order to approx- 
imate the spacetime between the numerical grid points. 
For larger times t, the geodesies behave as expected; in 
particular they are boosted to the right in order to evade 
the singularity and are attracted towards the left. For 
times sufficiently close to t = 0, however, the geodesies 
show an unexpected behavior; it seems that they are 
pushed away from the black hole region. 

Let us describe this phenomenon in more detail and 
relate it to certain components of the Riemann tensor, see 
Fig. |5J There we show a contour plot of the x-component 
of the vector field 

a" := R^ttx, 

defined from the Riemann curvature tensor. This quan- 
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Figure 6. Contour plot of a x as defined in the text. 

tity is the relative acceleration vector of two neighboring 
geodesies which both start in the i-direction and which 
are separated only in the x-direction; see [45[ . Its leading- 
order behavior for Gowdy spacetimes turns out to be 

a x = (1 - k 2 )r 2 + . . . 

This is compatible with the first plot in Fig. [5] and hence 
is the general behavior close to any Gowdy singularity 
with k < 1. It is, in particular, positive insight most of 
the black hole-like region in our case. In physical terms 
this means that such causal geodesies, which are directed 
away from the singularity, accelerate away from each 
other. Close to the intermediate point in Fig. H3 how- 
ever, we see that a x becomes negative. Here, such causal 
geodesies are rather accelerated towards each other when 
going away from the singular region. This is a phe- 
nomenon where the effective leading-order descriptions 
fails; in fact, the expansion of the Riemann curvature at 
the intermediate point is trivial since all terms of arbi- 
trary order vanish. However, the first of the two approx- 
imation schemes, which underlies the numerical compu- 
tations, is able to unveil this reliably. We see that the 
region where a x is negative coincides exactly with that 
part of the spacetimes where the geodesies in Fig. [5] show 
the previously described surprising behavior. This is cer- 
tainly a consequence of the particular choice of asymp- 
totic data, which "concatenate" a region with very strong 
gravitational fields to a region with vanishing gravita- 
tional fields, that the resulting gravitational fields have 
this surprising property. 

V. CONCLUDING REMARKS 

In this paper we have discussed a variety of ideas rel- 
evant to the Fuchsian method and the singular initial 
value problem for general second-order hyperbolic Fuch- 
sian equations and we have presented a perspective from 
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the physics standpoint. Our method allows us, on one 
hand, to check under which specific conditions a given 
leading-order description of the dynamics of evolution 
equations is consistent and, on the other hand, yields 
two kinds of approximation schemes which allow us to 
go beyond this effective description and, therefore, to 
approximate actual solutions with in principle arbitrary 
accuracy. Hence, we construct solutions with prescribed 
singular behavior; one of the methods yields expansions 
of arbitrary order and hence is useful for qualitative stud- 
ies; the other method approximates the actual solution 
by a sequence of solutions to the (standard) initial value 
problem and, hence, is particularly useful for numerical 
implementation. We stress that it is in general a challeng- 
ing problem to find the correct leading-order term for the 
singular initial value problem associated with any given 
evolution equation. We have demonstrated that the pro- 
posed canonical two-term expansion yields the correct 
leading-order behavior in our applications, at least after 
some additional considerations. 

Wc have applied this method to Gowdy spacetimes. It 
was discovered earlier that Gowdy solutions to Einstein's 
vacuum equations can be constructed by means of the 
above method. In the present paper, we have in particu- 
lar investigated the asymptotic dynamics of freely falling 
observers in such a spacetime. We confirm, as previously 
expected, that inhomogeneities do not effect the leading- 
order behavior of the geodesies; it is solely at a higher or- 
der that inhomogeneities play a significant role. In turns 
out that geodesies that are orthogonal to the symmetry 
orbits all have a limit point on the t — 0-surface. The 
same conclusions also hold for non-orthogonal geodesies 
if \k\ < 1. If k — ±1, however, the geodesies do not in 
general have a limit point, but can rather "swirl" towards 
the singularity. 

Let us mention here that, of course, it is not triv- 
ial to define geometrically what we mean by a "limit 
point on a curvature singularity" . Statements like "two 
geodesies approach the same point on the singularity" 
are in general rather coordinate dependent as well. For 
example, let us consider a family of timelike orthogonal 
geodesies on a Gowdy spacetime close to t — where k 
is smaller than one which emanate from the same spa- 
tial point x* at t — in the sense before; the singular 
initial value problem above indeed allows us to construct 
such geodesies. Then, locally for t > 0, we can intro- 
duce Gaussian coordinates where those geodesies repre- 
sent the timelines, and therefore, with respect to these 
new coordinates, each of these geodesies has a different 
limit point on the singularity. However, such a family 
Fi of geodesies is indeed geometrically distinct from any 
family F2 where the orthogonal geodesies approach dif- 
ferent spatial points at t = with respect to the original 
coordinates above. Let us consider the particle horizon 
in the sense of [if?] at any time > of an observer 
traveling away from the singularity at t = along one 
of the orthogonal timelike geodesies of Fi. Then it fol- 
lows that all the other geodesies of Fi intersect the past 



light-code of the observer at and are therefore never 
completely beyond the observer's particle horizon, irre- 
spective of how small i* is. The same situation leads to 
a different conclusion for the family F2; there is always a 
critical time such that any other given geodesic in F2 
does not intersect the past light-cone of the observer for 
all < T*. Hence, two future directed observers in the 
case of F2 cannot communicate with each other if they 
are too close to the singularity. But they always can in 
the case of Fi . 

This also has consequences for the quest for a rigor- 
ous formulation of the BKL conjecture. As discussed 
in [29], one imagines to introduce Gaussian coordinates 
locally close to the singularity. The claim is that the 
field equations effectively decouple along each coordinate 
timeline and each timeline evolves as an independent 
Mixmaster-like universe asymptotically. Now our discus- 
sion for the Gowdy spacetimes shows that in fact, each 
Gaussian coordinate timeline with respect to the family 
Fi approaches the same, in this case, Kasner universe 
asymptotically. The claim of the BKL conjecture there- 
fore only holds with respect to the family F2. Whether 
or not, it is possible to distinguish these two types of 
geodesic congruences a-priori, say, on the level of the ini- 
tial data for the initial value problem, is not clear in 
general for a system of equations as complicated as Ein- 
stein's field equations. 

Finally, we have studied a particular Gowdy solution 
numerically, that is, a solution with a black hole-like re- 
gion next to a Minkowski-like region exhibiting a Cauchy 
horizon. Of particular interest was the region where these 
two regions come together and interact. This led us to 
unexpected properties for the behavior of those causal 
geodesies which go into the intermediate point. It would 
be interesting to check whether similar phenomena are 
present for different types of asymptotic data. In short, 
we expect that the proposed theory and its generaliza- 
tions (e.g. [ill ]) will be useful for various problems arising 
in physics and, especially, general relativity. Applying 
the theory to more general classes of spacetimes is work 
in progress. 
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